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1. Introduction 

Systems which feature s-wave pairing between femrions of different species — say f and \. 
fermions — exhibit interesting phase structure as we vary the difference in their chemical poten- 
tials, <5/i = | (/if — Hi)/2\. For small enough <5/i, the free energy gain provided by Cooper pairing 
drives the system into a superfluid state. All fermions are paired up, and the number densities of 
t particles (w-j-) and J, particles are equal. As we increase 8fl, pairing is stressed since in the 
absence of pairing the corresponding Fermi surfaces would be split and the system would be po- 
larized, 8n = — ^ 0. At large enough 5/i the pairing is disrupted and the system exists in a 
normal phase. 1 

Stressed pairing of this kind shows up in a variety of physical systems. In quark matter at 
densities relevant to the phenomenology of neutron stars, the strange quark mass gives rise to a 
splitting between the effective quark chemical potentials [0]. In cold atoms, polarization can be 
created by trapping different numbers of fermions of the two hyperfine states that pair. Since the 
number of particles in the traps is fixed, polarization leads to phase separation []3|] which can be 
observed [Q] by looking at their density profile. 

In these proceedings, we will focus on the pairing between spin f and J, neutrons in the inner 
crust of neutron stars. Dipole interactions with the magnetic field split the chemical potentials of 
the two. At a critical value of the field, pairing breaks, and there is a transition to the normal phase 
from the superfluid phase. 

Calculating the critical field is relevant for the phenomenology of neutron-star crusts because 
the transport properties of the two phases are very different. The superfluid phase can transport heat 
efficiently via Goldstone bosons (although mixing between the superfluid Goldstone bosons and 
the lattice phonons suppresses this effect [|5[ ^]). Furthermore, the normal phase features gapless 
neutrons near the Fermi surface, and consequently has a large specific heat. Both these effects 
tend to increase the time scales [|7|] associated with heat transport in the inner crust if a significant 
fraction of the crustal neutrons are in the normal phase. 

The phase diagram as a function of 5/i is well understood in BCS theoiy. The normal phase is 
favored over the superfluid phase for 5/i = 8jX c > Aq/t/2 where Ao is the pairing gap for Sjl = 0. 
In a small window near 5jJ, c exotic phases called LOFF [j|] phases are theoretically favored. 

Similar results can be argued for "weak coupling" [Q: a regime where the separation between 
the particles is much greater than the scattering length. More formally, defining kp = (37i 2 (n^ + 
n^)) 1//3 , and the scattering length as a, the free energy of the normal and the superfluid phases can 
be written in terms of functions that become smaller as \kpa\ decreases, if kp a is small. 

For strongly coupled Fermi systems, however, we don't know of a controlled expansion tech- 
nique to calculate thermodynamic properties. Therefore, one needs to perform intensive numerical 
calculations to compute the energies of both the normal and the superfluid phase. 

In these proceedings we will present QMC calculations of the energy of the competing phases 
for three number densities n\ = 6.65 x 10~ 4 fm~ 3 , n% = 2.16 x 10~ 3 fm~ 3 , and 723 = 5.32 x 10~ 3 
fm -3 for different polarizations. For comparison, these correspond respectively to 0.004 2n sat , 
0.0135« sat , and 0.033?2 sa t where n sat = 0.16 fm~ 3 is the saturation density. The neutron-neutron 

'Here we are assuming that the masses of the two species are equal. A similar stress arises in systems of fermions 
with different masses (for eg. 
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x = N i /N T 

Figure 1: Ground-state energy per particle (in units of the free Fermi gas energy) scaled to the power of 
3/5 for normal spin-polarized neutron matter. Shown are QMC results at three different total densities 
rei,/i2,rt3 (as points: black circles, red squares, blue diamonds, respectively), along with fits to the Monte 
Carlo results using the set of functions of (x,p) given in Eq. 3.6. (as lines: black solid, red dashed, blue 
dotted, respectively). Also shown using hollow symbols are the results for an unpolarized superfiuid. 



interaction has a large scattering length (a ~ —18.6 fm) compared to the effective range r e , ~ 2.7 
fm. We see that at the densities we consider kpa is large in magnitude (kpa = 5.02, 7.44, and 
10.05 corresponding to the three densities which are all deep in the strongly coupled regime). For 
reference, we note that kp r e = 0.59, 0.88, and 1.19. 

Phase competition is analyzed very simply in the grand canonical ensemble. If we know the 
pressure of the competing phases as a function of the average chemical potential pL = + Hi)/ 2 
and the splitting 8fX, then the phase with the larger pressure at given }X wins. At the phase transition 
the pressure and /J. should match. (S/J. need not be the same between coexisting phases, as long as 
8/J. < A, where A is the superfiuid gap.) 

The direct evaluation of the pressure as a function of ji is very difficult because of the sign 
problem. For a fixed ru, n± a creative technique known as fixed node QMC [10] has been developed 
to calculate the energy of fermionic systems S(n^ Sn) accurately. We will describe this calculation 
next. 



2. Calculation of S 

The heart of the calculation is the accurate evaluation of the energy as a function of density 



and polarizations [ |11| , |12fl . We use the Diffusion Monte Carlo (DMC) method which is expected to 
be exact, modulo the fermion-sign problem. Schematically, this entails taking a variational wave 
function *¥y that encodes knowledge about the physics of the system (e.g. is it superfiuid, does 
it contain trimers, etc.). We first perform a Variational Monte Carlo simulation and then project 
out the exact, lowest-energy eigenstate *Po from the trial (variational) wave function *Fy. In other 
words, we treat the Schrodinger equation as a diffusion equation in imaginary time and evolve the 
variational wave function up to the limit of % — > °°. This can be seen by expanding *Py in terms of 
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the complete set *F/ of the exact eigenstates of the Hamiltonian (with eigenvalues £7): 

= ^e-^-^^o, lirnr -> 00 . 



(2.1) 



where Ej is called the trial energy. In practice, the DMC technique is implemented by discretizing 
X. 

As already implied, the trial wave function is critically important, in that it is used as the fixed- 
node constraint. In addition, it is commonly also used as an importance-sampling function. In 
this problem, we are interested in simulating non-superfluid systems with unequal populations of 
spin-up and spin-down neutrons. Thus, modulo a conceptually irrelevant (but numerically useful) 
Jastrow term, the wave function describes the particles as being in a free Fermi gas (two, to be 
precise, one for each spin-orientation). This formalism localizes all the correlations within the 
Jastrow functions, leaving the rest of the wave function as a product of two Slater determinants: 



<J> S (R) = ^[^ n {r\)^n{ri) ■ ■ ■ 0„(^)]^[0„(r 1 /)0„(r 2 ') . . . 0„ (/#/)] 



(2.2) 



where the unprimed (primed) indices correspond to spin-up (spin-down) neutrons, N^ + N'^ =N, 
srf is the antisymmetrizer, and the n (r^) are plane waves with periodic boundary conditions inside 
a box of volume L 3 . Using real numbers is often preferable numerically, so it is customary to use 
linear combinations of the eigenf unctions. Shown below is the result for 7 particles: 



1 

cos(2kx\ /L) 
cos(2ny\ /L) 
cos{2kz\ /L) 
sin(2jixi/L) 
sin(2ny\/L) 
sin(2nz\/L) 



1 

cos{2kx2/L) 
cos(2ny2/L) 
cos(2kzi/L) 
sin{2nx2 /L) 
sin(2ny2/L) 
sin(2Kzz/L) 



1 

cos(2kx-i jt) 
cos(2jtyj jt) 
cos(2nzi /L) 
sin(2nxy / L) 
sin{2ny-i /L) 
sin(2%zi /L) 



(2.3) 



The number 7 was chosen because it corresponds to a closed shell. Closed shells exhibit no ambi- 
guity in the selection of the momenta different particles are placed in. 

In our simulations of polarized normal neutron matter we used 33 + 0, 57 + 19, 57 + 33, 33 + 
27, and 33 + 33 particles. These correspond to relative fractions of x = 0, 0.333, 0.579, 0.818, and 
1, respectively. At each particle number combination, we performed simulations for three distinct 
total number densities [n = (N^+N^/L 3 ], as mentioned earlier. The results of these computations 
are given in Fig. [|, which shows the values of 



g{x,n) 



S, 



FG(t) 



3/5 



(2.4) 



where ^G(t) = f ^ 2m ' ( ?1 t) 5 ^ 3 ^ s tne ener gy of a non-interacting gas of a single species. 
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3. The free energy functions 

3.1 The unpolarized superfluid phase 

The energy and the gap for the unpolarized superfluid phase for n\, ti2 and 713 were calcu- 



lated in [13, h4p. A convenient parameterization for the data can be obtained by considering the 



dimensionless function gsF(n), defined by the relation 

o (t_2\2/3 

<?{n) = ^J—{ng SF {n)fl' , (3.1) 

where, $ is the energy density. The prefactors and the power of n in the definition is chosen so that 
for a system of non-interacting fermions, g = 1. Since g$ F is dimensionless, the density depen- 
dence in gsF can only arise in appropriate combinations with a and r e . In particular, it motivates 
considering the form 

gSF{n) = go,o + go,\r e n l/3 + g L0 —j— -\ (3.2) 

For the relevant densities, a good fit is obtained for go o = 0.942, go,ir e = 0.003MeV _1 , g\$/a = 
2.301MeV. For comparison, for the unitaiy Fermi gas {a — > °°, r e — > 0), g$F is independent of the 
density and is given by (2§) 3 / 5 , where % is the Bertsch parameter. For t, = 0.38 [15] calculated 
using similar techniques, (2£ ) 3 / 5 = 0.848, which is close to go.o- It is not surprising that go,o differs 
from the unitary value since this calculation is being done for neutron matter, and the interaction is 
more complicated. One expects deviations also because the largest kp r e is of the order of 1. 

With this parameterization, it is easy to calculate the pressure as a function of the chemical 
potential using standard thermodynamical relations. 

3.2 The normal phase 

The analysis of the normal phase is a little more involved because the energy depends both on 
the density n and the relative fraction x = n^/n^ where we take to be the majority species. A 

is 

3 (6tt 2 ) 2 / 3 



convenient generalization of Eq. 3. 1 



-(n t g(x,n)) 5 ^ , (3.3) 



5 2m 

which is equivalent to Eq. |2.4[ The chemical potential is: 

_d£ _ 3 2 / 3 7r 4 / 3 (^ T ) 2/3 g(x,n) 2 / 3 (^jf(x,n) + (l-x)jf(x,n)+g(x,»; 

dn Sn 2 4 / 3 m 
The chemical potential splitting is: 

dS 3^K^{^) 2/3 g(x,n)^((-x-l)^(x,n)+g(x,n)) 



5/i 



ddn 



2 4/3 



111 



(3.4) 



(3.5) 



We fit the values of g given in Fig. |] with a linear combination of the following functions, 

{( * +1) '^' (x+1)pl/3 ' (x + 3)2 ' 1 ^ 

(3.6) 

The x dependence is chosen so that at 8pL = for x = 1 . 
The pressure is simply P = —$(n, 8n) + fin + 8/J.Sn. 
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Figure 2: Pressure as a function of 8ji for fixed ji = 1.94 MeV. For djJ. > 8/i c = 0.53 MeV the unpaired 
phase (bold black curve) has a higher pressure. At 8jj. c the SF and the NP phases can coexist and for smaller 
8fi, the SF phase is favored. On the other extreme for 8fi = 4.18 MeV, x = 0.0005. We choose so that 
the density of the normal phase at 8jJ, c is «2- 



4. Phase competition 

For a given /», 8pL c is simply the point where the pressure curves of the unpolarized superfluid 
and the normal phases cross. For example, from the plot of the pressure as a function of SjJ. for 
ll = 1.94 MeV in Fig. g we conclude that 5/i c = 0.53 MeV. 

In Fig. |3] we show a plot of 8pL c /Ep as a function of }i, where Ep = k\j (2m) (m is the mass of 
the neutron and kp is the Fermi momentum corresponding to the density in the superfluid phase). 
The range of jj. is chosen such that the density of the superfluid varies from n\ at the smallest value 
to «3 at the largest. For reference we note that A/Ep increases gradually from 0.4 at n\ to 0.44 at 



«3 [13]. This can be seen as a competition between stronger pairing as \kpa\ increases weakened 
by increasing kpr e . The effect of the effective range on 8pL c /Ep is much more pronounced. We 
see that, not surprisingly, this makes it "easier" (in a dimensionless sense) to break Cooper pairing 
at higher density. The two curves give an estimate of the uncertainty associated with different 
interpolating functions of the same data. 



5. Conclusions 



To summarize, we have performed Diffusion Monte Carlo simulations for polarized neutron 
matter at low density. Low density implies that the interaction between neutrons is greatly sim- 
plified: three4}ody forces and spin-orbit & tensor interactions have a minimal effect, leading to a 
mainly s-wave interaction. We have included the s-wave channel of the high-quality phenomeno- 
logical Argonne vl8 potential, along with perturbative p-wave corrections. The computations we 
have performed were undertaken for normal neutrons of varying population numbers and are ex- 
pected to provide tight upper bounds to the true ground-state energies. 

Using these results we have extracted the critical chemical potential splitting 8pL c at which 
there is a transition from the superfluid state to the normal state as we increase 8fl. 
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H(MeV) 



Figure 3: Plot of 8fi/Ef as a function of ji. The two curves correspond to two choices of interpolating 
functions in x. The solid line (black) corresponds to the set of functions given in Eq. 3.6, up to cubic 
power in x. The dashed line (red) corresponds to a larger set including quartic powers. The vertical lines at 
increasing fi correspond to the chemical potentials at which the densities of the unpolarized superfluid phase 
are «i, «2 and «3 respectively. 



For a system with unequal number densities of f and 4- neutrons our calculations imply that 
the energetically favored state is that of phase separation into an unpolarized superfluid and a po- 
larized normal neutron gas. This situation directly parallels what is found in systems of ultracold 
fermionic atoms [[If]] which have no effective-range and p-wave contributions. The cold-atom sys- 
tem has also been probed experimentally and quantitatively well-described by QMC calculations, 
lending credence to the claim that our results for neutron matter are also quite accurate. Using our 
simulations and worked published earlier, we have extracted the critical relative fractions above 
which phase separation occurs (for details see [11]). For relative fractions above the critical values 
the competing phase that has been considered [ |12| ] is a homogeneous polarized superfluid state. 
The mixed phase is energetically favored compared to the polarized superfluid state, meaning that 
we have provided the lowest-energy configuration to date for the system under study. At relative 
fractions lower than the critical value, our analysis suggests that the system is entirely normal. 

The quenching of pairing, if it materializes in a neutron star sufficiently magnetized to polar- 
ize the neutrons, would have directly observable consequences, as superfluidity impacts the specific 
heat and more generally the thermal behavior of the star. This would also hold even for slightly 
smaller magnetic fields. More generally, our predictions are in principle also relevant to exper- 
iments with cold atoms which could in the future be performed for finite effective ranges in the 
laboratory. [|7|] Such experiments today make use of unequal spin populations as a matter of 
course. Importantly, our microscopic simulations can also impact the terrestrial phenomenology 
of nuclei. Skyrme-family energy-density functional theories of nuclei have traditionally been fit- 
ted to experimentally measured masses of nuclei, while also sometimes including infinite matter 
constraints. The latter in the past consisted of the behavior of nuclear matter, while more recently 
the neutron-matter pairing gap has also been used as an extra constraint. [18, [T^] Following this 
path to its logical consequences, a recently completed work [20] directly compares the neutron 
polaron binding energy (i.e. the binding of one spin-down impurity embedded in a sea of spin- 
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up neutrons) with the predictions following from a group of established Skyrme parametrizations. 
Extending this line of thinking, one could generalize such comparisons to finite relative fractions 
and thus provide a further constraint to nuclear density functionals: this all follows from the fact 
that our DMC results provide a dependable microscopic solution of the Schroedinger equation. 
Neutron-rich nuclei, in particular, are expected to be impacted by the physics of polarized neutron 
matter. 

There are many ways in which one could build on the calculations discussed in this contribu- 
tion. One could add an external periodic potential, which would stand for the periodic lattice of 
neutron-rich nuclei found in a neutron-star crust and could also teach us something about surface 
effects in neutron drops. The physics of the static response of polarized (or, for that matter, unpolar- 
ized) neutron matter is wholly uncharted from a microscopic perspective. Any predictions on this 
system could also be extended to the relevant setting of optical-lattice experiments, which could 
indirectly teach us something about the neutron case, even without an effective range. Similarly, 
one could extend calculations such as these to more species or even finite-temperature effects, in 
an effort to directly probe thermal effects in neutron-star crusts and cores. 
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